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We present a quantitative study of vorticity formation in peripheral ultrarelativistic heavy ion collisions at 
sfstm = 200 GeV by using the ECHO-QGP numerical code, implementing relativistic dissipative hydrodynam¬ 
ics in the causal Israel-Stewart framework in 3+1 dimensions with an initial Bjorken flow profile. We consider 
and discuss different definitions of vorticity which are relevant in relativistic hydrodynamics. After demonstrat¬ 
ing the excellent capabilities of our code, which proves to be able to reproduce Gubser flow up to 8 fm/c, we 
show that, with the initial conditions needed to reproduce the measured directed flow in peripheral collisions 
corresponding to an average impact parameter b = 11.6 fm and with the Bjorken flow profile for a viscous Quark 
Gluon Plasma with tj/s = 0.1 fixed, a vorticity of the order of some 1(T 2 c/fm can develop at freezeout. The 
ensuing polarization of A baryons does not exceed 1.4% at midrapidity. We show that the amount of developed 
directed flow is sensitive to both the initial angular momentum of the plasma and its viscosity. 


I. INTRODUCTION 

The hydrodynamical model has by now become a paradigm 
for the study of the QCD plasma formed in nuclear colli¬ 
sions at ultrarelativistic energies. There has been a consider¬ 
able advance in hydrodynamics modeling and calculations of 
these collisions over the last decade. Numerical simulations 
in 2+ID HI and in 3+1 D t2j-(7| including viscous corrections 
are becoming the new standard in this field and existing codes 
are also able to handle initial state fluctuations. 

An interesting issue is the possible formation of vorticity in 
peripheral collisions 18 • TT)J. Indeed, the presence of vortic¬ 
ity may provide information about the (mean) initial state of 
the hydrodynamic al evolution which cannot be achieved oth¬ 
erwise, and it is related to the onset of peculiar physics in the 
plasma at high temperature, such as the chiral vortical effect 
DU- Furthermore, it has been shown that vorticity gives rise 
to polarization of particles in the final state, so that e.g. A 
baryon polarization - if measurable - can be used to detect 
it I1 1 211131 . Finally, as we will show, numerical calculation 
of vorticity can be used to make stringent tests of numerical 
codes, as the T-vorticity (see sect. [II] for the definition) is ex¬ 
pected to vanish throughout under special initial conditions in 
the ideal case. 

Lately, vorticity has been the subject of investigations in 
refs. J9][l0l with peculiar initial conditions in cartesian coor¬ 
dinates, ideal fluid approximation and isochronous freezeout. 
Instead, in this work, we calculate different kinds of vortic¬ 
ity with our 3+1D ECHO-QGP [] code 0, including dissi¬ 
pative relativistic hydrodynamics in the Israel-Stewart formu¬ 
lation with Bjorken initial conditions for the flow (i.e. with 


1 The code is publicly available at the web site http://theory.fi.infn.it/echoqgp 


u x = u y = u n — 0), henceforth denoted as BIC. It should be 
pointed out from the very beginning that the purpose of this 
work is to make a general assessment of vorticity at top RHIC 
energy and not to provide a precision fit to all the available 
data. Therefore, our calculations do not take into account ef¬ 
fects such as viscous corrections to particle distribution at the 
freezeout and initial state fluctuations, that is we use smooth 
initial conditions obtained averaging over many events. 


A. Notations 

In this paper we use the natural units, with ti = c - K — 1. 
The Minkowskian metric tensor is diag(l, -1, -1, -1); for the 
Levi-Civita symbol we use the convention e 0123 = 1. 

We will use the relativistic notation with repeated indices as¬ 
sumed to be summed over, however contractions of indices 
will be sometimes denoted with dots, e.g. u ■ T ■ u = u^T^Uy. 
The covariant derivative is denoted as (hence d^g )lv = 0), 
the exterior derivative by d, whereas d t , is the ordinary deriva¬ 
tive. 


II. VORTICITIES IN RELATIVISTIC HYDRODYNAMICS 

Unlike in classical hydrodynamics, where vorticity is the 
curl of the velocity field v, several vorticities can be defined 
in relativistic hydrodynamics which can be useful in different 
applications (see also the review fl4]| ). 
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A. The kinematical vorticity 


This is defined as: 


(O^v — ^ (ffyMp d^Uy) — dyU p (£ u v ) 


(1) 


where u is the four-velocity field. This tensor includes both 
the acceleration A and the relativistic extension of the angu¬ 
lar velocity pseudo-vector oj fI in the usual decomposition of 
an antisymmetric tensor field into a polar and pseudo-vector 
fields: 


tUyjv — ^jlV/XT^d Lf ^ (Ajjlly A yllji) 

A m = 2 u)[i V u v = u'dyU p = Du p 

CO^ — — fyjpa-T td M 


( 2 ) 


where ep ypo - is the Levi-Civita symbol. Using of the transverse 
(to u) projector: 


A MV = / v - vdu\ 

and the usual definition of the orthogonal derivative 
Vp =: A ^d a — d^i — u^D, 

where D = u a d a , it is convenient to define also a transverse 
kinematical vorticity as: 

w pv = \ p A vo -af lr = ^(VyMp - VpMy) (3) 

Using the above definition in the decomposition <[2]» it can be 
shown that: 

Wpy = tfjvpo-of M°" (4) 

that is co A is the tensor formed with the angular velocity vector 
only. As we will show in the next subsection, only oj a shares 
the “conservation” property of the classical vorticity for an 
ideal barotropic fluid. 


B. The T-vorticity 


This is defined as: 

\ [dyC T«p) - <9pC Tu v )\ (5) 

and it is particularly useful for a relativistic uncharged fluid, 
such as the QCD plasma formed in nuclear collisions at very 
high energy. This is because from the basic thermodynamic 
relations when the temperature is the only independent ther¬ 
modynamic variable, the ideal relativistic equation of motion 
(s + p)Ajj = Vpp can be recast in the simple form (see e.g. 

113): 


id f^py — 


1 

2 


(TA V - VyT) = 0 


( 6 ) 


The above ([6| is also known as Carter-Lichnerowicz equa¬ 
tion m for an ideal uncharged fluid and it entails conserva¬ 
tion properties which do not hold for the kinematical vortic¬ 
ity. This can be better seen in the the language of differential 
forms, rewriting the definition of the T-vorticity as the exterior 
derivative of a the vector field (1-form) Tu, that is ft = d( 7 u). 
Indeed, the eq. ([6]i implies - through the Cartan identity - that 
the Lie derivative of Q along the vector field u vanishes, that 
is 


£ U Q. - u ■ dQ + d (u ■ D.) - 0 (7) 

because £1 is itself the external derivative of the vector field Tu 
and dd = 0. The eq. (|7]i states that the T-vorticity is conserved 
along the flow and, thus, if it vanishes at an initial time it will 
remain so at all times. This can be made more apparent by 
expanding the Lie derivative definition in components: 

U u fir = oo! lv - d lT idn ,rv - d a u v n ,ni = o (8) 

The above equation is in fact a differential equation for Q pre¬ 
cisely showing that if £2 = 0 at the initial time then £1 = 0. 
Thereby, the T-vorticity has the same property as the classical 
vorticity for an ideal barotropic fluid, such as the Kelvin cir¬ 
culation theorem, so the integral of O over a surface enclosed 
by a circuit comoving with the fluid will be a constant. 

One can write the relation between T-vorticity and kinemat¬ 
ical vorticity by expanding the definition ([5]): 

Opy = ^ [(SyT) Mp - (dpT) My] + T Wpy 

implying that the double-transverse projection of £1: 

Appear = £2py = Toj% 

Hence, the tensor oj a shares the same conservation properties 
of Q A , namely it vanishes at all times if it is vanishing at the 
initial time. Conversely, the mixed projection of the kinemat¬ 
ical vorticity: 


lfaiprtr = 

does not. It then follows that for an ideal uncharged fluid with 
w A = 0 at the initial time, the kinematical vorticity is simply: 

COflV — —( AfjUy — AyUfj) (9) 


C. The thermal vorticity 

This is defined as lEl: 

= ^(<3yy8p - <VL) (10) 

where ft is the temperature four-vector. This vector is de¬ 
fined as (1 /T)u once a four-velocity m, that is a hydrody¬ 
namic al frame, is introduced, but it can also be taken as a 
primordial quantity to define a velocity through u = (3/ yfjfi 
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m. The thermal vorticity features two important properties: 
it is adimensional in natural units (in cartesian coordinates) 
and it is the actual constant vorticity at the global equilib¬ 
rium with rotation HU for a relativistic system, where /? is 
a Killing vector field whose expression in Minkowski space- 
time is P/j - b M + m /Jv x v being b and m constant. In this case 
the magnitude of thermal vorticity is - with the natural con¬ 
stants restored - simply fujo/ksT where to is a constant angular 
velocity. In general, (replacing to with the classical vorticity 
defined as the curl of a proper velocity field) it can be read¬ 
ily realized that the adimensional thermal vorticity is a tiny 
number for most hydrodynamical systems, though it can be 
significant for the plasma formed in relativistic nuclear colli¬ 
sions. 

Furthermore, the thermal vorticity is responsible for the lo¬ 
cal polarization of particles in the fluid according to the for¬ 
mula lfl2l : 

n P (x,p) = -le upcrT (l - n F ) m p ‘ rI — (11) 

8 m 


if the transversely projected vorticity tensor co A initially van¬ 
ishes, so will the transverse projection Q A and <zr A and the 
kinematical and thermal vorticities will be given by the for¬ 
mulae <0 and ( [T4| respectively. Indeed, the T-vorticity O 
will vanish throughout because also its longitudinal projec¬ 
tion vanishes according to eq. (0. This is precisely what hap¬ 
pens for the usually assumed BIC for the flow at To, that is 
u x = u y - u n - 0, where one has to A = 0 at the beginning as 
it can be readily realized from the definition <[T]>. On the other 
hand, for a viscous uncharged fluid, transverse vorticities can 
develop even if they are zero at the beginning. 

It should be noted though, that even if the space-space 
components (x, y, rj indices) of the kinematical vorticity ten¬ 
sor vanish at the initial Bjorken time to, they can develop at 
later times even for an ideal fluid if the spatial parts of the 
acceleration and velocity fields are not parallel, according to 
eq. (|9]). The equation makes it clear that the onset of spatial 
components of the vorticity is indeed a relativistic effect as, 
with the proper dimensions, it goes like (a x v)/c 2 . 


which applies to spin 1/2 fermions, n F being the Fermi-Dirac- 
Juttner distribution function. 


n F 


1 

oPM-p-n/T + 2 


( 12 ) 


Similarly to the previous subsection, one can readily obtain 
the relation between T-vorticity and thermal vorticity: 


^ [ (d M T ^ u y ~ ( a > T ) u ^ Q, 


T 2 


(13) 


Again, the double transverse projection of m is proportional 
to the one of Q: 


A /jp A vo -tj7 ptr = m A v = = iw A 


whereas the mixed projection turns out to be, using eq. (13 1 


vPw„ A™ = -\v v T + — 

per 2t2 2t 


Again, for an ideal uncharged fluid with to A = 0 at the initial 
time, by using the equations of motion ([6]), one has the above 
projection is just A v /T and that the thermal vorticity is simply: 


V — y, (AfjUy - 


A v iiu) 


(14) 


A common feature of the kinematical and thermal vorticity 
is that their purely spatial components can be non-vanishing 
if the acceleration and velocity field are non-parallel, even 
though velocity is vanishing at the beginning. 


In the full longitudinally boost invariant Bjorken picture, 
that is u n = 0 throughout the fluid evolution, in the ideal case, 
as « A = 0, the only allowed components of the kinematical 
vorticity are co TX , to ry and to xy from the first eq. <0- The a> xy 
component, at r/ = 0, because of the reflection symmetry (see 
fig-0 in both the x and y axes, can be different from zero but it 
ought to change sign by moving clockwise (or counterclock¬ 
wise) to the neighbouring quadrant of the xy plane; for central 
collisions it simply vanishes. 

However, in the viscous case, more components of the vor¬ 
ticities can be non-vanishing. Furthermore, in more realistic 
3+1 D hydrodynamical calculations, a non-vanishing u n can 
develop because of the asymmetries of the initial energy den¬ 
sity in the x - rj and y - r\ planes at finite impact parameter. 
The asymmetry is essential to reproduce the observed directed 
flow coefficient t'i(y) in a 3+1D ideal hydrodynamic calcula¬ 
tion with BIC, as shown by Bozek ED, and gives the plasma 
a total angular momentum, as it will be discussed later on. 

In this work, we calculate the vorticities, and especially the 
thermal vorticity m by using basically the same parametriza- 
tion of the initial conditions in ref. Da. Those initial con¬ 
ditions are a modification of the usual BIC to take into ac¬ 
count that the plasma, in peripheral collisions, has a relatively 
large angular momentum (see Appendix A). They are a min¬ 
imal modifications of the BIC in that the initial flow velocity 
Bjorken components are still zero, but the energy density lon¬ 
gitudinal profile is changed and no longer symmetric by the 
reflection // —> -rj. They are summarized hereinafter. Given 
the usual thickness function expression: 


III. HIGH ENERGY NUCLEAR COLLISIONS 

In nuclear collisions at very large energy, the QCD plasma 
is an almost uncharged fluid. Therefore, according to pre¬ 
vious section’s arguments, in the ideal fluid approximation. 


X oo /-*oo 

d zn(x,y,z)= I dz- ”° - 

oo J -co 1 + e ( y x2+ y 2+ A-R)is 

(15) 

where ;?o = 0.1693 fm 3 , 6 — 0.535 fm and R = 6.38 fm 
are the nuclear density, the width and the radius of the nu¬ 
clear Fermi distribution respectively, the following functions 






and 
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Figure 1: Colliding nuclei and conventional cartesian 
reference frame. Also shown the initial angular momentum 
vector. 

are defined: 

r >= r -H 1 - ^f) 

where cr is the inelastic NN cross section, A the mass number 
of the colliding nuclei, and: 

T + (x t ) = T(x t + b/2) T.(x t ) = T(x t - b/2) (18) 

where x T = (x, y) is the vector of the transverse plane coor¬ 
dinates and b is the impact parameter vector, connecting the 
centers of the two nuclei. In our conventional cartesian refer¬ 
ence frame, the b vector is oriented along the positive x axis 
and the two nuclei have initial momentum along the z axis 
(whence the reaction plane is the xz plane) and their momenta 
are directed so as to make the initial total angular momentum 
oriented along the negative y axis (see fig. [TJ. The wounded 
nucleons weight function Wn is then defined: 

W N (x,y,rf) = 2 (T x (x,y)f-{ri) + T 2 (x,y)f+(j))) (19) 

where: 


(16) 

(17) 


f-(n) 


l 

-T} + T]„, 
7 TJ m 


V < -Tim 
-rim <ri<T] m 

n > r] m 


f+Oi) = 


0 

rj + rim 
2r) m 

1 


ri < -rim 
-rim <71 < rim 
ri > rim 


Finally, the initial proper energy density distribution is as¬ 
sumed to be: 


e(x,y, rj) = e 0 W(x,y, rj) H(rj), (20) 

where the total weight function W(x,y, rf) is defined as: 

(1 - a) W N (x,y, rj ) + an BC (x,y) 


W(x,y , rj) 


11 n ) W;v(0,0,0) + arise) 0,0' 


( 21 ) 


b=0 


and: 


//(//) = exp 


Tj / 


r) = \rj\ - T) f i at l2 ( 22 ) 


In the eq. (21 1 nsc(x,y) is the mean number of binary colli¬ 
sions: 


r>Bc(x,y) = cr m T+(x,y) T-{x,y) (23) 

and a is the collision hardness parameter, which can vary be¬ 
tween 0 and 1. 

This parametrization, and especially the chosen forms of 
the functions / ± , are certainly not unique as a given angular 
momentum can be imparted to the plasma in infinitely many 
ways. Nevertheless, as has been mentioned, it proved to re¬ 
produce correctly the directed flow in a 3+ID hydrodynami- 
cal calculation of peripheral Au-Au collisions at high energy 
ED, thus we took it as a good starting point. A variation of 
this initial condition will be briefly discussed in sect. |VII| Be¬ 
sides, the parametrization ( [20] ) essentially respects the causal¬ 
ity constraint that the plasma cannot extend beyond // = ybeam- 
Indeed, at V*nn = 200 GeV ybeam — 5.36 while the 3 cr point 
in the gaussian profile in eq. (j22]) lies at r/ = r]t\ M /2+3cr rl - 4.4. 

The free parameters have been chosen following ref. ED, 
where they were adjusted to reproduce the data in Au-Au col¬ 
lisions at V%n = 200 GeV. They are reported in table [l] 

We have run the ECHO-QGP code in both the ideal and 
viscous modes with the parameters reported in table [I] and the 
equation of state reported in ref. j20l . The impact parameter 
value b — 11.57 was chosen as, in the optical Glauber model, 
it corresponds to the mean value of the 40-80% centrality class 
(9.49 < b < 13.42 fm M2T1 ) used by the STAR experiment for 
the directed flow measurement in ref. 123. The initial flow 
velocities u x , u y , u ,] were set to zero, according to BIC. The 
freezeout hypersurface - isothermal at Tf„ - 130 MeV - is 
determined with the methods described in refs. 0(231. 


IV. QUALIFICATION OF THE ECHO-QGP CODE 

To show that our code is well suited to model the evolution 
of the matter produced in heavy-ion collisions and hence to 
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Parameter 

Value 

V^WV 

200 GeV 

a 

0. 

e 0 

30GeV/fm 3 

v in 

40 mb 

T 0 

0.6 fm /c 

T\flat 

1 

<T„ 

1.3 

Tfo 

130 Me V 

b 

11.57fm 

77,,, ideal 

3.36 

T] m viscous 

2.0 

r//s 

0.1 


Table I: Parameters defining the initial configuration of the 
fluid in Bjorken coordinates. The last two parameter values 
have been fixed for the last physical run. 


carry out our study on the development of vorticity in such an 
environment, we have performed two calculations, referring 
to an ideal and viscous scenario respectively, providing a very 
stringent numerical test. 

Before describing these tests, it should be pointed out that 
the vorticities components are to be calculated in Bjorken co¬ 
ordinates, whose metric tensor is g MV = diag(l, -1,-1, -r 2 ), 
hence they do not all have the same dimension nor they are 
adimensional as it is desirable (except the thermal vorticity, as 
it has been emphasized in Sect. [EI]i. For a proper comparison 
it is better to use the orthonormal basis, which involves a fac¬ 
tor r when the ?/ components are considered. Moreover, the 
cumulative contribution of all components is well described 
by the invariant modulus, which, for a generic antisymmetric 
tensor A MV is: 

A 2 = = 2 [A 2 xy - A 2 x - A 2 t + (A 2 x +A 2 y -A 2 T )/r 1 ]. (24) 


Furthermore, we have always rescaled the T-vorticity by 1/T 2 
in order to have an adimensional number. Since the T-vorticity 
has always been determined at the isothermal freezeout, in 
order to get its actual magnitude, one just needs to multiply it 

b y T %' 


of the grid resolution (the boundaries in x, y, tj being fixed) 0 
As it is expected, the normalized T-vorticity decreases as the 
resolution improves. 

Because of the relation ( fl3] i, the residual value at our best 
spatial resolution of 0.15 fm can be taken as a numerical error 
for later calculations of the thermal vorticity. 
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Figure 2: (color online) Mean of the absolute value of 
T-vorticity components, divided by T 2 , at the freeze-out as a 
function of the grid resolution. 
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A. T-vorticity for an ideal fluid 

Since the fluid is assumed to be uncharged and the initial T- 
vorticity £1 is vanishing with the BIC, it should be vanishing 
throughout, according to the discussion in sect. 0- However, 
the discretization of the hydrodynamical equations entails a 
numerical error, thus the smallness of Q in an ideal run is a 
gauge of the quality of the computing method. In fig. [2] we 
show the mean of the absolute values of the six independent 
Bjorken components at the freeze-out hypersurface, of the T- 
vorticity divided by T 2 to make it adimensional, as a function 


Figure 3: (color online) Contour plot of fl r ,/r7’ 2 at the 
freeze-out hypersurface at y = 0. 


2 It should be pointed out that, throughout this work, by mean values of 

the vorticities we mean simple averages of the (possibly rescaled by 1/r) 
Bjorken components over the freezeout hypersurface without geometrical 
cell weighting. Therefore, the plotted mean values have no physical mean¬ 
ing and they should be taken as descriptive numbers which are related to 
the global features of vorticity components at the freeze-out. 
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B. Gubser flow 


A very useful test for the validation of a numerical code 
of relativistic dissipative hydrodynamics is the explicit solu¬ 
tion of Israel-Stewart theory of a Bjorken flow with an az- 
imuthally symmetric radial expansion EH3, the so-called 
Gubser flow. Indeed, this solution provides a highly non¬ 
trivial theoretical benchmark. 

For the sake of clarity, we briefly summarize the main steps 
leading to the analytical solution, to be compared with the 
numerical computation. In the case of a conformal fluid, 
with p = e/3 EOS, the invariance for scale transformations 
sets the terms entering the second-order viscous hydrody¬ 
namic equations. The additional requests of azimuthal and 
longitudinal-boost invariance, constrain the solution of the 
hydrodynamic equations, which has to be invariant under 
SO(3 ) 9 ® SO(l, 1) ® Z 2 transformations. To start with, one 
defines a modified space-time metric as follows (with usual 
Bjorken coordinates, 77 being the spacetime rapidity): 


ds 1 = T 2 ds 2 = t 2 




? 

= T 




which can be viewed as a rescaling of the metric tensor: 


ds 2 > dS 2 = ds 2 /! 2 <^> gp V > gp V = gfiv/T 2 . 

It can be shown that ds 2 is the invariant spacetime interval of 
dSj ® R, where dS 3 is the three-dimensional de Sitter space 
and R refers to the rapidity coordinate. It is then convenient 
to perform a coordinate transformation (q is an arbitrary pa¬ 
rameter setting an energy scale for the solution once one goes 
back to physical dimensionful coordinates) 


sinhp = - 


1 - g-(r- - r-) 
2 qT 


tan# = 


2 qr 


1 + q 2 (r 2 - r 2 ) ’ 


(25) 


after which the rescaled spacetime element ds 2 reads 

ds 2 = dp 2 - cosfrp (dd 2 + sini 9d(/) 2 ) - drf. (26) 


most general form actually admits further terms that were de¬ 
rived for a system of massless particles in refs. USED), valid 
for the case of a conformal fluid with s - 3p: 



In the case of the Gubser flow in Eq. (28 1 , due to the traceless 
and transverse conditions ft^ = 0 and iqjv v = 0 , one has simply 
to solve the two equations (n' 1 ’ 1 = 7 r w /sT) 


1 dT 2 1 __ u 

—-1— tanh p = —n" tanh p 

T dp 3 3 


(29) 


and ( 77 /s = 77 /s, being the ratio dimensionless) 


tr 


fGfa’) 2 . a„h P 
dp 3 


+ 7r w = —%■ tanh p. 
3 sT 


(30) 


The solution can be then mapped back to Minkowski space 
through the formulae: 


T = f/T, 


Kpv — 


1 dx a dx@ „ 
r 2 dx^ dx v7Ial3 


(31) 


In fig. [4] we show the comparison between the Gubser analyt¬ 
ical solution and our numerical computation for the temper¬ 
ature T and the components n xx , n xy and 7r w of the viscous 
stress tensor respectively, at different times. The initial energy 
density profile was taken from the exact Gubser solution at the 
time r = 1 fm/c. The simulation has been performed with a 
grid of 0.025 fm in space and 0.001 fm in time. The shear vis¬ 
cosity to entropy density ratio was set to 77 /s = 0.2, while the 
shear relaxation time is t r = 5q/(s + p). The energy scale is 
set to q = 1 far 1 . As it can be seen, the agreement is excellent 
up to late times. 


C. T-vorticity for a viscous fluid 


The full symmetry of the problem is now manifest. SO(l, 1) 
and Z 2 refer to the usual invariance for longitudinal boosts 
and 77 —> -77 inversion, while SO(3) ? reflects the spherical 
symmetry of the rescaled metric tensor in the new coordinates. 
In Gubser coordinates the fluid is at rest: 

Up = 1, uq — U(f) — Ujj — 0. (27) 

The corresponding flow in Minkowski space can be obtained 
taking into account both the rescaling of the metric and the 
change of coordinates 

<9F„ 

U/J — T— ~Uy , 

dxf* 

where x 4 ‘ = (p, 6, (p, 77 ) and F = (r, r, c/> , 77 ). Other quantities 
such as the temperature or the viscous tensors require the so¬ 
lution of the following set of hydrodynamic equations (their 


Unlike for an ideal uncharged fluid, T-vorticity can be gen¬ 
erated in a viscous uncharged fluid even if it is initially van¬ 
ishing. Thus, the T-vorticity can be used as a tool to estimate 
the numerical viscosity of the code in the ideal mode by ex¬ 
trapolating the viscous runs. 

A comment is in order here. In general, in addition to stan¬ 
dard truncation errors due to finite-difference interpolations, 
all shock-capturing upwind schemes are known to introduce 
numerical approximations that behave roughly as a dissipative 
effects, especially in the simplified solution to the Riemann 
problems at cell interfaces If30ll . It is therefore important to 
check whether the code is not introducing, for a given reso¬ 
lution, numerical errors which are larger than the effects in¬ 
duced by the physics. We refer to the global numerical errors 
generically as numerical viscosity. 

We have thus calculated the T-vorticity for different physi¬ 
cal viscosities (in fact 77 /s ratios), in order to provide an upper 
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Figure 4: (color online) Comparison between the semi-analytic solution of the Gubser viscous flow (solid lines) and the 

numerical ECHO-QGP computation (dots). 


bound for the numerical viscosity of ECHO-QGP in the ideal 
mode. The mean value of the T-vorticity is shown in fig. [5] 
and its extrapolation to zero occurs when \tj/s\ < 0.002 which 
is a very satisfactory value, comparable with the one obtained 
in ref. ©. The good performance is due to the use of high- 
order reconstruction methods that are able to compensate for 
the highly diffusive two-wave Riemann solver employed 13. 


V. DIRECTED FLOW, ANGULAR MOMENTUM AND 
THERMAL VORTICITY 

With the initial conditions reported at the end of the Sect. m 
we have calculated the directed flow of pions (both charged 
states) at the freezeout and compared it with the STAR data 
for charged particles collected in the centrality interval 40- 
80% Il22l . Directed flow is an important observable for sev¬ 
eral reasons. Recently, it has been studied at lower energy 
ED with a hybrid fluid- transport model (see also ref. Il32l ). 
At V^nn = 200 GeV, it has been calculated with an ideal 
3+ID hydro code first by Bozek m. Herein, we extend the 
calculation to the viscous regime. 

The amount of generated directed flow at the freezeout de¬ 


pends of course on the initial conditions, particularly on the 
parameter rj m (see Sect. as shown in fig. [ 6 ] The directed 
flow also depends on rj/s as shown in fig.[7]and could then be 
used to measure the viscosity of the QCD plasma along with 
other azimuthal anisotropy coefficients. It should be pointed 
out that, apparently, the directed flow can be reproduced by 
our hydrodynamical calculation only for -3 < y < 3. The de¬ 
pendence of V! (y) on 77 ,„ and 77 /s makes it possible to adjust the 
77 ,,, parameter for a given 77 /s value. This adjustment cannot be 
properly called a precision fit because, as we have mentioned 
in the Introduction, several effects in the comparison between 
data and calculations have been deliberately neglected in this 
work. However, since our aim was to obtain a somewhat re¬ 
alistic evaluation of the vorticities, we have chosen the value 
of 77 ,,, for which we obtain the best agreement between our 
calculated pion vi(y) and the measured for charged particles 
in the central rapidity region. For the fixed 77 /s = 0.1 (ap¬ 
proximately twice the conjectured universal lower bound) the 
corresponding best value of rj m turns out to be 2.0 (see fig. [ 8 ]». 

It is worth discussing more in detail an interesting rela¬ 
tionship between the value of the parameter 77 ,,, and that of 
a conserved physical quantity, the angular momentum of the 
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T|/S 


Figure 5: (color online) Mean of the absolute values of 
Q /n ,/7' 2 components at the freeze-out hypersurface as a 
function of tj/s. Note that the 0 A7; , O v); , Q n; have been 
multiplied by 1/r. Upper panel: log scale. Lower panel: 
magnification of the region around zero viscosity. 


plasma, which, for BIC is given by the integral (see Appendix 
A for the derivation): 


P 


= - T °f 


dv dy d/f x s(x , y, rj ) sinh rj 


(32) 


Figure 6: (color online) Directed flow of pions for different 
values of r\ m parameter with ///,v = 0.1 compared with STAR 
data f22l . 



Y (rapidity) 


Figure 7: (color online) Directed flow of pions for different 
values of tj/s with ij m — 2.0 compared with STAR data |[22l . 


sharp spheres. In our conventional reference frame, the initial 
angular momentum of the nuclear overlap region is directed 
along the y axis with negative value and can be written as: 


Since r\ m controls the asymmetry of the energy density distri¬ 
bution in the rj - x plane, one expects that J y will vary as a 
function of rj, n . Indeed, if the energy density profile is sym¬ 
metric in rj , the integral in eq. ( [32} vanishes. Yet, for any finite 
rj m 4 0, the profile (20 1 is not symmetric and J y 4 0 (looking 
at the definition of /+ and f it can be realized that only in the 
limit rj m —> oo the energy density profile becomes symmetric). 
The dependence of the angular momentum on rj m with all the 
initial parameters kept fixed is shown in fig. [9] For the value 
rj m = 2.0 it turns out to be around 3.18 x 10 3 in h units. 

It is also interesting to estimate an upper bound on the an¬ 
gular momentum of the plasma by evaluating the angular mo¬ 
mentum of the overlap region of the two colliding nuclei. This 
can be done by trying to extend the simple formula for two 


I 


J' = I d.rdy w(x,y)(T + - T )x- 


/snn 


(33) 


where T± are the thickness functions like in eq. 18 and 
min(n(x + b/2 , y, 0), n(x - b/2, y, 0)) 


w(x,y ) = 


max(n(x + b/2,y, 0), n(x - b/2,y, 0)) 


is the function which extends the simple product of two 6 
functions used for the overlap of two sharp spheres. Note that 
the oj(x, y) is 1 for full overlap (b-0 ) and implies a vanish¬ 
ing angular momentum for very large b (see fig.flO} (see also 
ref. El). 

At b = 11.57 fm the above angular momentum is about 
3.58 x 10 3 in h units. This means that, with the current 
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Figure 8: Directed flow of pions at rj/s = 0.1 and i) m = 2.0 
compared with STAR data l22l . 



Figure 9: Angular momentum (in h units) of the plasma with 
Bjorken initial conditions as a function of the parameter rj m . 



0 2 4 6 8 10 12 14 

b (fm) 

Figure 10: (color online) Estimated angular momentum (in Ti 
units) of the overlap region of the two colliding nuclei (solid 
line) and total angular momentum of the plasma according to 
the parametrization of the initial conditions (dashed line), as 
a function of the impact parameter. 
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parametrization of the initial conditions, for that impact pa¬ 
rameter about 89% of the angular momentum is retained by 
the hydrodynamical plasma while the rest is possibly taken 
away by the corona particles. 

With the final set of parameters, we have calculated the 
thermal vorticity m. As it has been mentioned in Sect. [II] 
this vorticity is adimensional in cartesian coordinates) and it 
is constant at global thermodynamical equilibrium ED. e.g. 
for a globally rotating fluid with a rigid velocity field. In rel¬ 
ativistic nuclear collisions we are far from such a situation, 
nevertheless some thermal vorticity can be generated, both in 
the ideal and viscous case. This is shown in figs. m and [12] 
It can be seen that the generated amount of thermal vor¬ 
ticity has some non-trivial dependence on the viscosity. Par¬ 
ticularly, as it is apparent from fig. 12 the m xr] component - 
which is directed along the initial angular momentum - has 
a non-vanishing mean value whose magnitude significantly 
increases with increasing viscosity. Its map at the freeze- 


Figure 11: (color online) Mean of the absolute value of 
thermal vorticity covariant components at the freeze-out as a 
function of rj/s. Note that the m xrj , vj yn , m Tr] have been 
multiplied by 1 /r. 


out, for a fixed value of the y coordinate y = 0, is shown 
in fig. 13 where it can be seen that it attains a top (negative) 
value of about 0.05 corresponding to a kinematical vorticity, 
at the freezeout temperature of 130 MeV, of about 0.033 c/fm 


= 10 


i22 c -l 


In this respect, the Quark Gluon Plasma would 


be the fluid with the highest vorticity ever made in a terres¬ 
trial laboratory. However, the mean value of this component 
at the same value of rj/s = 0.1 is of the order of 5.4 x 10~ 3 , 
that is about ten times less than its peak value, as shown in 
fig. 12 This mean thermal vorticity is the consistently lower 
than the one estimated in ref. ED (about 0.05) with the model 
described in refs. (9] jTO] implying an initial non-vanishing 
transverse kinematical and thermal vorticity m A . This reflects 
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Figure 12: (color online) Mean values of thermal vorticity 
components at the freeze-out as a function of rj/s. Note that 
the m xn , Tu yri , m Tn have been multiplied by 1 /r. 
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Figure 13: (color online) Contour plot of 1 /r-scaled qx 
covariant component of the thermal vorticity, ct, ;a /t over the 
freeze-out hypersurface for y = 0 , q/s=0. 1 , q m =2.0. 


in a quite low value of the polarization of A baryons, as it will 
be shown in the next section. 


interesting feature of this relation is that it makes it possible 
to obtain an indirect measurement of the mean thermal vor¬ 
ticity at the freezeout by measuring the polarization of some 
hadron. For instance, the polarization of A baryons, as it is 
well known, can be determined with the analysis of the angu¬ 
lar distribution of its decay products, because of parity viola¬ 
tion. The polarization pattern depends on the momentum of 
the decaying particle, as it is apparent from eq. (34 1 . 

The formula ( |34| makes sense only if the components of 
the integrand are Minkowskian, as an integrated vector field 
yields a vector only if the tangent spaces are the same at each 
point. Before summing over the freezeout hypersurface we 
have then transformed the components of the thermal vortic¬ 
ity from Bjorken coordinates to Minkowskian by using the 
known rules. The thus obtained polarization vector W(p) is 
the one in the collision frame. However, the polarization vec¬ 
tor which is measurable is the one in the decaying particle rest 
frame which can be obtained by means of the Lorentz trans¬ 
formations: 


n 


o 

o 



m 


p n 

m 


no = n- 


p n 

e(m + e) 


P 


(35) 


In figure[T4]we show the A polarization vector components, 
as well as its modulus, as a function of the transverse momen¬ 
tum P 7 - for p z = 0 expected under the assumptions of local 
thermodynamical equilibrium for the spin degrees of freedom 
maintained till kinetic freezeout. It can be seen that the po¬ 
larization vector has quite an assorted pattern, with an over¬ 
all magnitude (see fig. [14} panel (a)) hardly exceeding 1% at 
momenta around 4 GeV. As expected, the y component is pre¬ 
dominantly negative, oriented along the initial angular mo¬ 
mentum vector and a magnitude of the order of 0.1%. Indeed, 
the main contribution to the polarization stems from the longi¬ 
tudinal component 11 ^, with a maximum and minumum along 
the bisector \p x \ = \p y \. 

The obtained polarization values are - as expected - consis¬ 
tently smaller than those estimated in ref. fl3ll (of the order 
of several percent with a top value of 8-9%) with the already 
mentioned initial conditions used in refs. O M. This is a 
consequence of the much lower value of the implied thermal 
vorticity, as discussed in the previous section. Also, the iH 
pattern is remarkably different, with different location of max¬ 
ima and minima. 


VI. POLARIZATION 


As it has been mentioned in the Introduction, vorticity can 
result in the polarization of particles in the final state. The 
relation between the polarization vector of a spin 1 /2 parti¬ 
cle and thermal vorticity in a relativistic fluid was derived in 
ref. tH and reads: 


m P ) = 


1 Jk dI,pV( 1 - n F ) p^fr^dyfip 
8m fz d£,ip A n F 


(34) 


where n F is the Fermi-Dirac-Juttner distribution function ( fj~2l > 
and the integration is over the freeze-out hypersurface X. The 


VII. CONCLUSIONS, DISCUSSION AND OUTLOOK 

To summarize, we have calculated the vorticities developed 
in peripheral (b — 11.6 fm) nuclear collisions at V^nn = 200 
GeV (b - 11.6 fm) with the most commonly used initial con¬ 
ditions in the Bjorken hydrodynamical scheme, by using the 
code ECHO-QGP implementing second-order, causal, rela¬ 
tivistic dissipative hydrodynamics. An extensive testing of the 
high accuracy and very low numerical diffusion properties of 
the code has been carried out, followed by long-time simula¬ 
tions (up to r = 8 fm/c) of the so-called viscous Gubser flow, a 
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Figure 14: (color online) Magnitude (panel a) and components (panels b,c,d) of the polarization vector of the A hyperon in its 

rest frame. 


stringent test of numerical implementations of Israel-Stewart 
theory in Bjorken coordinates. 

We have found that the magnitude of the 1 /t x - ij com¬ 
ponent of the thermal vorticity at freezeout can be as large as 
5 x 10 2 and yet its mean value is not large enough to produce 
a polarization of A hyperons much larger than 1%, which is a 
consistently lower estimate in comparison with other recent 
calculations based on different initial conditions. We have 
found that the magnitude of directed flow, at this energy, has 
an interestingly sizeable dependence on both the shear viscos¬ 
ity and the longitudinal energy density profile asymmetry pa¬ 
rameter //„, which in turn governs the amount of initial angular 
momentum retained by the plasma. 

The fact that in 3+ID the plasma needs to have an initial an¬ 
gular momentum in order to reproduce the observed directed 
flow raises the question whether the Bjorken initial condition 
u n — 0 is a compelling one or, instead, the same angular mo¬ 
mentum can be obtained with a non trivial u' 1 and with a suit¬ 
able change of the energy density profile. For a testing pur¬ 


pose, we have run ECHO-QGP with an initial profile: 

u' 1 = - tanhAx sinh(y b eam - H) (36) 

T 


which meets the causality constraint (see Appendix B). It is 
found that the directed flow is very sensitive to an initial it' 1 . 
For a small positive value of the parameter A = 5x 10 4 fm 1 
corresponding to a 7 V = 3.32 x 10 3 , keeping all other parame¬ 
ters fixed, the directed flow exhibits two slight wiggles around 
midrapidity (see fig. 15 i which are not seen in the data. For 
a very small negative value of the parameter A — -5 x 10 4 
fm' 1 , corresponding to J y = 3.08 x 10 3 , the directed flow in¬ 
creases while approximately keeping the same shape as for 
A — 0 around midrapidity. However, more detailed studies 
are needed to determine whether a non-vanishing initial flow 
velocity is compatible with the experimental observables. 


We plan to extend this kind of calculation to different cen¬ 
tralities, different energies and with initial state fluctuations in 
order to determine the possibly best conditions for vorticity 
formation in relativistic nuclear collisions. 
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APPENDIX A - ANGULAR MOMENTUM 

The calculation of the total angular momentum of the 
plasma can be done provided that initial conditions are such 
that energy density falls off rapidly at large \rj\. This condi¬ 
tion, which is met by the profile in eq. ( [20] ), indeed implies 
that a boundary exists where the angular momentum density 
tensor (that is the integrand below) vanishes and the following 
integral is conserved: 

F v = J dS,; ( x /J T Av - x v T‘ ifJ ) (37) 

where 2 is any spacelike hypersurface extending over the re¬ 
gion where the angular momentum density vanishes. The ob- 
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vious choice for Bjorken-type initial conditions is the hyper¬ 
surface T = Tq. 

It should be stressed that a vector (or tensor) integral is 
meaningful in flat spacetime only if the components are the 
cartesian ones. Hence, for the hypersurface r = to, the in¬ 
tegration variables are conveniently chosen to be the Bjorken 
ones, but the components of the stress-energy tensor as well 
as the x vector will be cartesian. Since the only non vanishing 
component of the angular momentum in our conventional ref¬ 
erence frame is P, orthogonal to the reaction plane, we can 
write: 


P = / 31 = J d£,i [. x 3 T a1 - x'T n \. (38) 

Finding the hypersurface measure dXj in cartesian compo¬ 
nents, but expressed through Bjorken variable, requires some 
reasoning. First, one has to remind that: 


dX, ( = dX n A (39) 

where n is the unit vector normal to the hypersurface r = To 
which is readily found to be (cartesian covariant components): 

nn = (cosh tj, 0,0, - sinh q) (40) 

Now, since Bjorken coordinates are time-orthogonal (g T! = 0) 
and with g TT = 1, the invariant spacetime measure dQ can be 
factorized into the product of the infinitesimal “time” dr and 
the infinitesimal measure of the orthogonal hypersurface dX: 

dQ = dxdX 


At the same time: 

dQ = Vi&idT d.v dy dq = t dr dx dy d// 


whence: 


dX = t dx dy dq (41) 


Using eqs. (39 1 , (40 1 and ( |41 | i, eq. ( [38] ) can be written ; 
P —T J dx dy dr/ [cosh^(x 3 r 01 - x^ 03 ) 


• sinh/ 7 (x 3 r 31 - x 


*r 33 )] 


(42) 


At the time t = to, the stress-energy tensor is supposedly the 
ideal one and there is no transverse velocity, so that 7’ 01 = 
T 31 = 0, while T 33 - (s + p)u z u z + p and T 03 = (s + p)u°u z . 
Pluggin these expressions into the ( |42| along with the trans¬ 
formation equation: 


t-rcoshr/ x-x y-y z = Tsinh/7 (43) 


one finally gets: 


P 


where: 


=to J 'dx dy d^ x cosh ij(e + p)u°u z 

+ sinh i] [(e + p)u z u z + /?)]] 


(44) 


m° = cosh rj u T + t sinh r/ u n 

u z — sinh rj u T + t cosh rj iP (45) 


being u T — Vl + Pu'i 2 . In the case of Bjorken initial condi¬ 
tions with uP = 0 and u T = 1, the eq. (|44| boils down to: 


P 


= - T °f 


dx dy d;; e(x, y, q) x sinh q 


(46) 


APPENDIX B - CAUSALITY CONSTRAINTS 


The inequality expressing the causality constraint in the hy- 
drodynamical picture of relativistic heavy ion collisions is that 
the initial longitudinal flow velocity must not exceed the ve¬ 
locity of beam protons v z < t'beam (assuming vanishing initial 
transverse velocity): 





1 + Vz \ 

1 -vJ 


| log(w° + U z )\ < y b eam 


(47) 


By using the transformation rules (45 i: 


log(w° + it) = log [(cosh 77 + sinh/7)(i< T + rw' 7 )] 

= log je^ Vl + t 2 uP + Tu n )J 
= q + log( Vl + PuP + Tu n ) 

= q + asinh( tm 7 ) < y be am ( 48 ) 


the inequality ( [47} becomes: 

I q + msinh(Tu n )\ < y bea m 


which can be solved for u n : 


- sinh(y beam + q)<u rl <- sinh(y beam - q) (49) 

T T 


The form (361 fulfills the above inequality. 











